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We examine Gaussian-basis Monte Carlo method (GBMC) introduced by Corney and Drum- 
mond. This method is based on an expansion of the density-matrix operator p by means of 
the coherent Gaussian-type operator basis A and does not suffer from the minus sign problem. 
The original method, however, often fails in reproducing the true ground state and causes sys- 
tematic errors of calculated physical quantities because the samples are often trapped in some 
metastable or symmetry broken states. To overcome this difficulty, we combine the quantum- 
number projection scheme proposed by Assaad, Werner, Corboz, Gull and Troyer in conjunction 
with the importance sampling of the original GBMC method. This improvement allows us to 
carry out the importance sampling in the quantum-number-projected phase-space. Some com- 
parisons with the previous quantum- number projection scheme indicate that, in our method, 
the convergence with the ground state is accelerated, which makes it possible to extend the 
applicability and widen the range of tractable parameters in the GBMC method. The present 
scheme offers an efficient practical way of computation for strongly correlated electron systems 
beyond the range of system sizes, interaction strengths and lattice structures tractable by other 
computational methods such as the quantum Monte Carlo method. 

KEYWORDS: Monte Carlo method, strongly correlated electron systems, Hubbard model, quantum- 
number projection 



1. Introduction 

Ground state properties of strongly correlated elec- 
tron systems are challenging subjects in condensed mat- 
ter physics. From the numerical point of view, there ex- 
ist many numerical algorithms, such as the exact diag- 
onalization method, the auxiliary-field quantum Monte 
Carlo (AFQMC) method, ^ the density matrix renor- 
malization group (DMRG) method 5 ^ and the path- 
integral renormalization group (PIRG) method. 6 - 11 ) Al- 
though the exact diagonalization of the Hamiltonian ma- 
trix gives accurate results, the tractable system size is 
severely limited. The AFQMC method can treat larger 
systems and has been applied to various correlated sys- 
tems. In some systems such as doped Mott insulators 
and the Mott insulators with geometrical frustration ef- 
fects, however, the AFQMC method often suffers from 
the negative sign problem which causes large statisti- 
cal errors by the cancellation of positive and negative 
Monte Carlo samples. The DMRG method offers prac- 
tically exact results without suffering from the negative 
sign problem. However, tractable lattice systems are re- 
stricted to one-dimensional configurations owing to the 
spatial renormalization process. The PIRG method is 
a powerful sign-free numerical technique for correlated 
electron systems and has been applied to various sys- 
tems beyond the tractable range of the above numer- 
ical methods. A main practical limitation in the PIRG 
method comes from the extrapolation procedure to reach 
the results for the full Hilbert space. The truncation error 
depends on the system size as well as on the interaction 
strength. 

Gaussian-basis Monte Carlo (GBMC) method has 



been proposed as an alternative quantum Monte Carlo 
method which does not involve any sign problem. 12, 13 ) 
This method is based on a representation of the density- 
matrix operator p by making use of the non-Hcrmitian 
Gaussian-type operator basis A. The Gaussian represen- 
tation is a natural generalization of the positive- P phase- 
space method which is often used in the area of quantum 
optics. 14 ' 15 ) As well as classical phase-space variables 
like (x,p), the Gaussian-basis representation utilizes the 

one-particle "Green's function" nij = Tr 



Tr 



= Tr 



and the stochastic weight 



fl as the phase-space variables, where cj (c$) is a Fermion 
creation (annihilation) operator of the i-th mode. In this 
method, we solve a Fokker-Planck equation with respect 
to the phase-space variables A = (fl,n,m,m + ), which 
is constructed by a mapping from an operator Liouville 
equation of p. One of the phase-space variables works 
as a weight of the importance sampling and for any 
two-body Hamiltonian, Q remains positive definite. Thus 
there exists no explicit manifestation of the negative sign 
problem. However, in many parameter regions, especially 
in the low-temperature region, the numerical results ob- 
tained by the GBMC method often show systematic er- 
rors. 16 ' 17 - 1 It will turn out below that this deviation is 
concerned with the "spontaneous symmetry breaking " . 

Assaad et at have used a quantum-number projec- 
tion scheme to overcome the deviations in the low- 
temperature region. 16, 17 ) They have proposed to project 
the density matrix onto given quantum numbers of the 
ground state after the sampling is completed and have 
reproduced accurate ground states in some parameter 
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regions. In their method, however, the convergence with 
the ground state becomes slower with the increase of the 
interaction strength, which determines the practical lim- 
itations. 

In this study, we combine the quantum-number projec- 
tion scheme concurrently with the importance sampling 
of the original GBMC method to reflect the amount of 
the overlap with the projected sector in the sampling 
weight. This allows us to perform the importance sam- 
pling with respect to the projected distribution. The cru- 
cial point is that the efficiency of the importance sam- 
pling is improved because the sampling weight reflects 
not only the energy but also the overlap with the pro- 
jected state, i.e., the overlap with the state which re- 
tains the same quantum numbers with the ground state. 
Thanks to the efficient sampling, the tractable parame- 
ter region becomes wider than the previous method re- 
ported by Assaad et al. 16 ' 17 * 1 Moreover, our method al- 
lows us to analyze the projected distribution directly. By 
using this advantage, the relation between the numerical 
convergence and the behavior of the projected distribu- 
tion is also reported. We show benchmark analysis up to 
10 x 10 lattices on the square lattice as well as up to the 
relative interaction strength U/t = 15 for the Hubbard 
model, which indicate the applicability and efficiency of 
the present method. 

The organization of this paper is as follows. Section 
2 gives an introduction of the GBMC method for the 
general Fermion systems by following the formulation in 
Refs. 12, 13, which we supply for the self-contained de- 
scription. In §3, we explain implementations of the Monte 
Carlo procedure for the Hubbard model. The improve- 
ments of the GBMC method by quantum-number pro- 
jections are shown in §4. In the last part of §4, we dis- 
cuss the practical limitation in the applicability of the 
GBMC method. Section 5 is devoted to summary and 
discussions. 

2. Gaussian-Basis Monte Carlo Method 

Gaussian-basis Monte Carlo (GBMC) method is a nu- 
merical method which makes use of the mapping between 
operator equations of motion and stochastic evolution 
equations of generalized phase-space. 12 ' 13 ^ In order to 
make an exact mapping, we introduce a complete set 
of Gaussian-type operators A, which is typically non- 
Hermitian. This basis set allows us to expand any phys- 
ical density-matrix operator p(r) in terms of the phase- 
space variables X as 



p(r) = / dXP(X, r)A(A), 



(1) 



where r is real or imaginary time, P(X, r) is the expan- 
sion coefficient, and dX is the integration measure. Here, 
-P(A) t) can always be chosen positive as in Appendix B 
and hence can be regarded as a probability distribution. 
In the GBMC method, it is the distribution P(X, r) that 
is sampled stochastically. 



2.1 General Gaussian basis 
2.1.1 Notation 

Before defining the Gaussian basis, we summarize 
the notation which will be used. Consider a M-mode 
Fermionic system characterized by the creation and an- 
nihilation operators and Cj , with anticommutation re- 
lations 

[di, c]]+ = 6 itj , [c tl Cj\ + = [cj, c]]+ = 0, (2) 

where i,j = l,..., M. We define a M-mode column vec- 
tor of the annihilation operators and its Hcrmitian con- 
jugate row vector as c and c\ respectively. In order to 
define the general Gaussian basis in a compact form, we 
introduce an extended-vector notation 



c 



ct = (ct,c T ). 



(3) 



Throughout the paper, we use the bold-type notation for 
M-mode vectors or matrices and the underline notation 
for 2M-mode extended vectors or matrices. 

For products of operators, we define a normal and an 
antinormal ordering operators denoted by : : and { }, 
respectively. The normal ordering operator : : reorders 
so that all the creation operators are put to the left of 
the annihilation operators, e.g., : :— —c^Ci. Simi- 
larly, the antinormal ordering operator { } reorders so 
that all the annihilation operators are put to the left 
of the creation operators, e.g., {cjci} — — c,cj. More 
generally, in the case of a nested product, the outer 
ordering operator does not reorder the inner one, e.g., 



{: c k Cj : Ci} — di : Cfec] := -C;c]cfc. The sign changes are 
necessary because of the anticommuting nature of the 
Fermion operators. 

2.1.2 Definition of the Gaussian basis 

Using the above notation, a general Gaussian operator 
is defined as : 



A (A) = SlPi\a] : exp - a' 1 /2)c 



where J is an extended unit matrix: 

—I 
I 



a is an extended covariance: 



n 



m 1 



m 
I n 



and the vector parameter A is defined as: 
A = (ft, n, m, m + ). 



(4) 



(5) 



(6) 



(7) 



Here, n is a M x M matrix which corresponds to nor- 
mal Green's function, while m and m + are two indepen- 
dent M x M antisymmetric matrices which correspond 
to anomalous Green's functions. Pfaffian of the antisym- 
metrized covariance g_ A appears so as to satisfy the nor- 
malization condition 



Tr[A(A)] = n. 



(8) 



Here, q_ A is constructed by moving each row in the lower 
half rows every after the row with the same indices in 
the upper half, and moving each column in the right half 
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columns every before the column with the same indices 
in the left half, i.e., 



du 



an 

Cll 

cm i 



biM 
diM 



CllM 
ClM 



bMM O.MM 
d-MM cmm 



(9) 



The general Gaussian operator A(A) itself may corre- 
spond to a density-matrix operator under the conditions 
that ■n) = n, m+ = and the eigenvalues of the ma- 
trix n lie in the interval [0, 1] . However, we do not restrict 
the Gaussian operator to such conditions and it allows 
the Gaussian basis to be an overcomplcte set which can 
expand any physical density-matrix operator. 

2.1.3 Properties of Gaussian basis 

From the above definition, the general Gaussian oper- 
ators satisfy some important identities. 12,13 ^ First, the 
Gaussian operator operated by the Fermion operator 
c and $ can be associated with differentiations of the 
Gaussian operators with respect to their parameters: 



5 f A : 



s <9A 



{c:c f A :} 



-a A 



< n 9A 



(10) 



(11) 



(12) 



Second, the traces of ladder operators with respect to 
the Gaussian operator can be analytically taken by using 
Grassmann coherent states: 



Tr 



Tr 



didjA 



Sin,- 



Qrrii 



Tr 



44 A 



(13) 
(14) 



To construct a mapping, we first substitute the expan- 
sion in Eq. (1) into the Liouville equation (15) to get 

/ dP< dt ^ ^ d - = J P( -' *)^&)^- ( 19 ) 

Second, using the differential properties in Eqs. (10-12), 
one can transform the superoperator ZA(A) into a differ- 
ential operator £A(A). We next apply partial integration 
to get, provided that boundary terms vanish, 

dP(X, t) 



J A(A) dP ^' ^ d ± = J A{X)jC'P(X,t)dX, 



(20) 



where C is reordered form of C. Note that from Eq. (10- 
11) £' contains derivatives only up to the second order 
for any two-body Hamiltonian. As a sufficient solution 
for Eq. (20) , a Fokker-Planck equation of Ito type is ob- 
tained: 

jP{\,t) = C'P(\,t) 



^E^^A) 



P{h,t). (21) 



The imaginary-time evolution equation of the density- 
matrix operator then boils down to the Fokker-Planck 
equation, which is in practice solved by integrating nu- 
merically the corresponding stochastic differential equa- 
tions (SDE). 

3. Gaussian Representation for Hubbard Model 

3. 1 Mapping 

We consider the following Hubbard Hamiltonian. 



H = 



N 



N 



(22) 



2.2 Time evolution 

The time evolution of a density-matrix operator is de 
termined by the Liouville equation 

dp(t) 



dt 



Lp(t), 



(15) 



where L is a Liouville superoperator. For a real time 
evolution, the superoperator is given by the commutator 
with the Hamiltonian: 



Lp(t) 



H,p(t) 



(16) 



In the case of calculating the equilibrium state at r — 
1/ksT, the imaginary time evolution of a density-matrix 
operator is determined by the equation 



dp{r) 

dT 



(17) 



Therefore the superoperator is given by the anticommu- 
tator with the Hamiltonian: 

Lp(r) = ~\H,p(r)\_ (18) 



-p for i = j, 

—t for being a nearest-neighbor pair, 
otherwise, 

(23) 

where i and j represent the lattice points, c\ a {ci a ) the 
creation (annihilation) operator of an electron with spin 
a on the i-th site, tij the transfer integral between the 
i-th site and the j-th site, U the on-site Coulomb in- 
teraction, p the chemical potential and N the number 
of the lattice sites. Although we treat only this simplest 
Hubbard model, the formation can easily be extended to 
a more general form including transfers for further-site 
pairs and/or intersite Coulomb interactions. 

Although Eq. (22) is a standard representation of the 
Hubbard Hamiltonian, it is necessary that the sign of the 
interaction term is negative so that the diffusion matrix 
.D(A) in Eq. (21) being positive definite. 12 ' 13 ) Thus, we 
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(24) 



2JV 



— T X yfl X y + \5 ar) 



N 



,n (iri),{ii7') 



x,y=l 



a. a' 



- n (ir,),{i V ') n (i<T).Xicr')\ Vr/ri'Vo 



(25) 



l xy 



where 

coordinates of site and spin, i.e., x 
operators and cj are defined as 



c* x c y and suffices x and y denote both the 

(i,a). The vector 



^1 — (^Tp ^ij.) 



(26) 
(27) 



The 2 AT x 2N extended hopping matrix T is defined as 
Tij = Ti + Nj + N = tij and Tfj+jv = ^i+JVj = 0, where 
i, j = 1, ■ ■ ■ , N. The matrix cr z denotes the 2 component 
of the Pauli matrix. Since the Hubbard model conserves 
the total particle number, we use the number-conserving 
subset of the general Gaussian operator to expand the 
density-matrix operator: 

-c+ [2/+(n T -iT 1 ]£ 



A(0, n) = Q dct(i - ro) : e 



(28) 



where n is a 2N x 2 AT matrix. 

The Gaussian operator consists of an overcomplctc set 
and it can expand any physical density-matrix opera- 
tor with positive coefficients. In the following sections, 
we express the parameters of the Gaussian operator as 
A = (ft, n). Similarly to the case of the general Gaussian 
operators in Eqs. (10-12), the number-conserving Gaus- 
sian satisfies the differential identities: 



A 



'■xy 1 



An, 



dk 

9n zw 



,jk -\- Ti XW {J) Z y Tlzy) 



dk 

dn z1 , 



(29) 



(30) 



(31) 



The trace of the Gaussian operator itself is Tr[A] = SI 
and the trace of any ladder operators can be calculated 
by Wick's theorem. For instance, we obtain 

Tr[A4cy = fln x 



(32) 

Ti'IAcJ.CycJyC^:] = [n xy n wz + n xz (5 W y — n wy )] . (33) 

To obtain the ground state of the system, one may con- 
sider the imaginary-time evolution of the density-matrix 
operator 

dp l j 

In the GBMC method, instead of solving the Liouville 
equation above, one solves generalized Langevin equa- 
tions by making use of the mapping between the Li- 
ouville equation and the stochastic equations. To this 



t = l/k B T. 



(34) 



end, we expand the density-matrix operator as p — 
J dAP(A, r)A(A). Then the Liouville equation becomes 



J d\k(\) dP{ ^ T) = J d\P(\,T) j-i [#,A(A)] + j . 

+ (35) 



The differential identities of the Gaussian operator en- 



able us to transform — | 
form: 



H, A(A) 



into a differential 



if, A 



x,y 

2 ^ J \ y 

■i x.v \ 



_d 



xy 

d 2 



i x,y 
w.z 



dflxydflujz 



where 



if(n) =Tr[A(n)J3]/Tr[A(n)], (37) 
A = ^n(T -UM)(I -n) 

+ i(i -n)(T -UM)n, (38) 



B (i) 

xy 



- \j^^Z a aa' n x,{ia'){S(i a ),y ~ "(iff),!,)) 



(39) 



C xy = \j~-^^2 a aa>( S x,(ia>) ~ "x,(iff') )"(iff) ,y 



(40) 



M, 



i/Vr),(jff') - % ^^^(^((T^Cr^, - (7^,(7^,). 
rj,rj' 



(41) 

Partial integration, under the assumption that boundary 
terms vanish, yields the Fokker-Planck equation for the 
probability distribution P(X, r): 



dP(X,r) 
dr 



x,y 



o 2 



xy 



#(») tfW 

U I l xv LS I v W z 



+ 



_ c (i) C (i) 



r)n 8n xy wz 



P(A,r). (42) 



In the actual calculation, instead of solving this equation 
directly, we solve the Ito-type Langevin equations with 
respect to the parameters of the Fokker-Planck equation 
which reproduce the distribution of P(X, t): 18 ) 

dtl = -QH(n)dT, (43) 
dn = -Adr + ^B«dWi +J2 cil)dW i ( 44 ) 
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where dW and dW are Wiener increments which sat- 
isfy (dWi) = (dWD = {dW l dW' j ) = and (dW t dW J } = 
(dW(dW^) = 5 tl dr . 

Any expectation values of physical observables are 
evaluated by using the trace properties in Eqs.(32) and 
(33). Let O be a general observable consisting of the lad- 
der operators, then the expectation value of O becomes 



(O) 



Tr[/3Q] 
Tr[p] 



JdXP(X, r)Tr 


A(A)0 


/ d\P(\, r)Tr 


A(A)" 



/ dXP{X, r)fiO(n) 



(45) 



Jd\P(X,T)Q 

In the GBMC method, the integration with the weight 
P(X, t) is achieved alternatively by summing up over all 
the walkers of the Langevin equations (43) and (44), i.e., 

E. 



(0) = 



E. n< 



(46) 



We now regard £1 as the weight of the importance sam- 
pling in the Monte Carlo procedure. Note that from 
Eq. (43), the formal solution of the weight £1 becomes 



£1(t) = exp 



f dT'H{n{ T ')) 
Jo 



(47) 



Since "Green's function" n(r) and H{n) are always real, 
the weight £1 remains positive. Hence the negative sign 
problem does not appear. 

3.2 Numerical integration 

When integrating the Langevin equations, one has to 
be careful about the type of the SDEs. Since Eq. (44) 
is Ito-type SDE, the numerical integration must be done 
by Ito integration. 18 ) Here, we introduce two schemes 
of the numerical integration. The simplest one is Euler- 
Maruyama scheme: 20 ^ 



x i+1 = x t + A(xi)Ar + B( Xi )AWi 



(48) 



This scheme is faster than any other scheme but is not 
stable in general. For a more stable integration, we use 
a semi-implicit iterative scheme 19 ^ : 

x l+1 =x l + A(x l+1 )AT + B(x l )AW i . (49) 

To solve the SDE, we make a first guess Xi + \ by Euler- 
Maruyama scheme (48). Then, is substituted into 
the drift term of (49) iteratively until a self-consistent 
solution is found. For the parameter values of the Hub- 
bard model we have chosen the time step At = 0.001, 
then only a few iterations are needed because the initial 
guess from the Euler-Maruyama scheme is already close 
to the final solution. 

3.3 Sampling method 

For an efficient calculation, the importance sampling 
is needed. From Eq. (46), Q can be regarded as a weight. 
Thus we can construct an importance sampling method 
with respect to 0. Corncy and Drummond use the 
branching method for importance sampling. 12 ' 13 > 21 ) The 
branching method works by cloning the samples whose 
weights are large and by killing whose weights are small. 



Assaad et al. also use a similar reconfiguration method 
but their method keeps total population constant. 16 ' 22 ) 
Here we propose another method which we call "suc- 
cessive Metropolis method" (see Fig. 1). In contrast 



T2 



1, n ^1 fj Q j se c h arl g e is accepted by 




r < 



9 Q 1 



Q 



warming up 



Steady Markov Chain 



randomly sample from this chain 



O O- 



o- 



o o 



i> 1 



r < 



o 



Q 



warming up 



O 



T 

Fig. 1. sampling method 

with the usual Metropolis method, this method allows 
us to evolve r successively. After a certain number of 
time steps NAt, Monte Carlo samples are stored by the 
Metropolis algorithm with the following conditions. 

• choose one starting point (Q^(ti), n^(ri)) 

and Wiener increments dWi 

■ evolve NAt steps to get (Qj(r2), nj(r2)) 

• choose another starting point (ti),^+i^ (ti)) 

- ( tr y) 

and another Wiener increments dW i+1 

■ evolve NAt steps to get (nf^ (t 2 ) , {t 2 )) 

■ select the new sample with the condition that 
(tt i+1 (T 2 ),n i+1 (T2),dWi +1 ) 

'(<T ) (r 2 ),n^ ) (r 2 ),d<T ) ) 

11 ' ^ n 4 (7^)/n;(Ti) 
k (O i (r 2 ) , rii (r 2 ) , dWi ) otherwise, 

where r is a uniform random number distributed in [0, 1). 
After a sufficient number of warming-up steps, the stored 
samples constitute a steady Markov chain which can be 
regarded as the new starting points of further time evo- 
lutions. 
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3.4 Systematic deviation 

Here, we demonstrate some elementary results in the 
case of the two-site Hubbard model at U/t = 4 and n = 1 
under the open boundary condition. Figure 2 shows the 
total energy, the specific heat C v and the charge suscep- 
tibihty X c = jk((N 2 ) - (N) 2 ) , where TV = £ i<T 4c i<7 . 
All the numerical results show excellent agreement with 
the exact diagonalization result. 




1.5 



0.9 



0.6 



0.3 
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0.12 
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'I' 1 

A (b)C V 


■ i 




Exact 




GBMC -B- 


1 2 


3 4 5 


T 





(c) X c - 






Exact 




GBMC -0- 



2 3 
T 



Fig. 2. (a) Energy of 2 X 1 lattice with U/t = 4 and n = 1 as a 
function of inverse temperature t. Solid curve is obtained from 
the exact diagonalization. (b) Specific heat as a function of tem- 
perature T obtained from numerical derivative of (a) . (c) Charge 
susceptibility as a function of temperature T. 



However, simulation results deviate from the exact di- 
agonalization results if the lattice size or the strength 
of the on-site interaction U becomes extremely larger. 



Here, as an example we demonstrate the results for the 
case of the two-site Hubbard model at U/t = 100 and 
n = 1. Figure 3 shows the total energy and the to- 



0.01 

oh 

-0.01 



<s 2 > 



-0.02 
-0.03 
-0.04 
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GBMC : 



Ei - 
E □ 



<S 2 > A 







10 
% 



15 



20 



1.4 
1.3 

1.2 
1.1 
1.0 
0.9 
0.8 
0.7 



Fig 



3. Energy and total spin on 2 X 1 lattice at U/t = 100 and 
n = 1 as functions of inverse temperature t. The squares and 
the triangles represent the GBMC results of the total energy 
and the total spin, respectively. Solid and dashed lines represent 
the exact value of the energy at the ground state Eq and the first 
excited state Ei, respectively. 



tal spin (S 2 ) for U/t = 100 and n = 1, where the to- 
tal spin operator S = (S x , S y , S z ) is defined as S a — 
I J2f =1 c\a a Ci, (a — x, y, z). As is seen from Fig. 3, the 
energy obtained by the GBMC method is located just 
at the middle between the ground state and the triplet 
first excited states. This means that the GBMC method 
reproduces the state which is represented by the super- 
position of the ground state and the triplet first excited 
states. Indeed, the expectation value of the total spin is 
(S 2 ) = 1, which is the middle point between the singlet 
state and the triplet states. Here, the ground state of 
the two-site Hubbard model at n = 1 is known to be 
represented as 



'1 + 



1 




(|^2> + |& 3 » - (50) 



where r = At/U and 



1 



^ sl ) = _(|U)-|U)) 



(51) 



|^2> = |U,0>, |<M = |o,U>. (52) 

Note that for r <C 1 , the ground state can be represented 

as 

l*o> *|&i> = ^(1 U> - I U»- (53) 

On the other hand, the triplet first excited states are 
represented as 



i^i> = ^(ira> + u,t» 



>I2 



) = IT,T), |&3> = |U>. 



(54) 
(55) 
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The GBMC result of S z is nearly zero in the whole range 
of t (not shown), and the GBMC method converges with 



|0i 8 > ± |&t> oc | U> or | |,T>. 



(56) 



This means that samples obtained by the GBMC method 
are trapped in a symmetry broken states as | f, J.) or | J. 
, t), because the electron hopping is prohibited owing to 
the energy loss by the strong Coulomb repulsion. 

When the system size becomes larger, the systematic 
deviation caused by the same reason as the case of the 
two-site Hubbard model at U/t — 100 occurs at a rela- 
tively small U/t. Here, for example we demonstrate the 
results for the case of the Hubbard model on the 4x1 
ring with n = 1 under the periodic boundary condition 
in the x direction. Hereafter, the 4x1 lattice results are 
all obtained from the same boundary condition. Figure 
4 shows the total energy in the case of U/t = 1 and 
U/t = 4. In the case of U/t = 1, the total energy agrees 
with the exact diagonalization result, while in the case 
of U/t — 4, the total energy deviates from the result of 
the exact diagonalization systematically. As is shown in 
the Table I, 4 x 1 lattice Hubbard model with U/t = 4 
and n = 1 has triplet exicited states with the energy 
E = -1.806424. We conclude that the GBMC samples 
are trapped in a symmetry broken state constructed from 
a linear combination of the singlet ground state and the 
triplet excited states. To confirm this, we have calculated 




Fig. 4. (a) Energy of 4 X 1 lattice with U/t = 1 and n = 1 as a 
function of inverse temperature r. Solid curve is obtained from 
the exact diagonalization. (b) Energy of 4x 1 lattice with U/t = 4 
and n = 1. 



Table I. Energy spectra of 4 X 1 lattice with U/t = 4 and n = 1 
under the periodic boundary condition in the onc-dimcnsional 
direction obtained by the exact diagonalization. 



Energy 


Number of Degeneracy 


-2.102748 


1 


-1.806424 


3 


-1.068140 


1 



the expectation value of the total spin {S 2 ). As is shown 
in Fig. 5, the total spin has a nonzero value and these 
overlaps with excited S > sectors cause the systematic 
deviation. 



<s 2 > 




Fig. 5. Expectation value of total spin as a function of inverse 
temperature r for the same case as Fig. 4 (b). 



3.5 Power-law tails 

To investigate the reason of the systematic deviation 
in detail, we calculate the distributions of the parameters 
of the Fokker-Planck equation (42) to analyze whether 
or not the boundary term in the partial integration of 
Eq. (35) exists. 

First, we calculate the distribution of the weight ft. As 
is seen form Fig. 6, the Monte Carlo step dependence of 
P(f2) implies that the distributions of both U/t = 1 and 
U/t = 4 have upper limits of SI, even though they have 
broad peak structures. To analyze this cutoff of P(f2) 
in detail, we also calculate the integrated distribution 
defined by 



Q(fi) = 1 - / P(fi')dn'. 
Jo 



(57) 



As we see from Fig. 7, the integrated distributions of both 
U/t = 1 and U/t = 4 have cutoffs around SI ~ 10 53 and 
SI ~ 10 124 , respectively. Thus there exists no boundary 
term with respect to the weight Q. 

Next, we calculate the distribution of Green's function 
n. Figure 8 shows the distribution of Green's function. 
The abscissa N denotes the square root of the sum of 
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Fig. 6. (color online): Distribution of weight f2 at r = 20 for the 
4x1 lattice at n = 1 in case of (a) U/t = 1 and (b) = 4. 
Red (solid), green (dashed) and blue (dotted) curves represent 
the distribution obtained by 1.28 X 10 6 , 2.56 X 10 6 and 3.84 X 10 6 
Monte Carlo steps, respectively. 
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Fig. 7. (color online): Integrated distribution of weight Q at r = 
20 for the case of 4 X 1 lattice at n = 1 in case of (a) U/t = 1 and 
(b) U/t = 4. Red (light), green (dashed) and blue (dark) curves 
represent the distribution obtained by 1.28 X 10 6 , 2.56 X 10 6 and 
3.84 X 10 6 Monte Carlo steps, respectively. 



each squared Green's function element, 



N = 




(if-), (jcr)' 



(58) 



and the ordinate P(N) denotes its probability distri- 
bution. Similarly to the case of P(£l), the distributions 
P(N) of both U/t — 1 and U/t = 4 have the upper limits. 
Thus there is no boundary term with respect to Green's 
function n, either. However, as is seen from Fig. 8, the 
distribution tails of both U/t = 1 and U/t = 4 show 
power-law-like behaviors, i.e., P(N) oc N~ r , below their 
cutoffs. If the exponent of the power law r is less or equal 
to p + 1, the p-th moment of Green's function diverges, 
i.e., 

r N p P(N)dN oc J N p -' r dN -► oo , if r < p + 1. 

(59) 

To estimate the exponent of the power law in detail, we 
have also calculated the integrated distribution defined 

as 

Q(N) = 1-1 P(N')dN'. (60) 
Jo 

Figure 9 shows the integrated distribution of both U/t = 
1 and U/t = 4. From the tails of Q(N), we obtain the 



power-law exponent of P(N), r = 5.42±0.03 for U/t = 1 
and r = 3.84±0.01 for U/t = 4, respectively. To make the 
energy be well defined, the power-law exponent r must 
be larger than three, since the energy is the second-order 
moment of Green's function (see Eq. (59)). In the case 
of U/t = 4, the exponent r is larger than three, which 
also supports the absence of the boundary terms in the 
partial integration of Eq. (35). From the analysis of the 
distribution P(f2) and P(N), we conclude that the sys- 
tematic deviations observed in the results of the orig- 
inal GBMC method are caused not from the boundary 
terms but from the trap to the quasi-stable states. In the 
next section, we introduce the quantum-number projec- 
tion method which can remove the systematic deviation. 



4. Quantum-Number Projection 

Generally speaking, quantum many-body systems 
have several symmetries inherent in a Hamiltonian such 
as translational symmetry, U(l) symmetry, SU(2) sym- 
metry, point group symmetry of the lattice, etc. Although 
these symmetries are sometimes broken in the thermody- 
namic limit, they must be preserved in finite size systems. 
In actual numerical calculations, however, these symme- 
tries are not always preserved in restricted Hilbert space 
and the numerical calculation often suffers from system- 
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Fig. 8. (color online): Distribution of Green's function at r = 20 
for the 4x1 lattice at n = 1 in case of (a) U/t = 1 and (b) 
U/t = 4. Red (solid), green (dashed) and blue (dotted) curves 
represent the distribution obtained by 1.28 X 10 6 , 2.56 X 10 6 and 
3.84 X 10 6 Monte Carlo steps, respectively. 



atic errors. 

One of the most promising device to restore these sym- 
metries is the quantum- number projection which has 
been used successfully in the framework of the path- 
integral renormalization group method. n ) Also in the 
framework of the GBMC method, Assaad et al. have used 
the quantum-number projection. 16 ' 17 ) They proposed to 
project the density matrix onto given quantum numbers 
of the ground state after the ordinary GBMC sampling 
is performed and reproduced accurate ground states in 
some parameter regions. In this paper, we call their 
method post-projected GBMC (GBMC-PS) method. 

In this section, we first review the mathematical frame- 
work of the quantum-number projection method, then 
introduce an alternative method for performing the 
quantum-number projection. In this method, we com- 
bine the quantum-number projection scheme in conjunc- 
tion with the importance sampling of the original GBMC 
method. This allows us to perform the importance sam- 
pling with respect to the quantum-number-projected dis- 
tribution, which makes it possible to treat wider param- 
eter region than the previous studies. 16, 17 - ) In this paper, 
we call this method pre-projected GBMC (PR-GBMC) 
method. 



10 u 
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Fig. 9. Integrated distributions at r = 20 for the 4x1 lattice 
at n = 1. Dashed curve represents the distribution of U/t = 1, 
whereas solid curve is that of U/t = 4. The logarithmic fit to 
the tails leads to the exponent r = 5.42 ± 0.03 for U/t = 1 
and r = 3.84 ± 0.01 for U/t = 4. Both data are obtained from 
3.84 X 10 6 Monte Carlo steps. 



4-1 Unitary transformation of a Gaussian operator 

Before introducing the quantum-number projectors, 
we define a unitary transformation of a Gaussian opera- 
tor which will be used in all the projectors. For any Her- 

mitian matrix h(x), a Gaussian operator is transformed 
agi6.iT) 



Ac 1 h(x)c 



A(fi,n) = A(n(x),n(x)), 



where 



fi(as) = n det \{e ih{x) - I)n T + / 



/ - n(x) = \ \(e ih ^ - I)n T + I 



(61) 



(62) 



(I-n). (63) 



To prove Eqs.(62) and (63), we first introduce several 
identities of the Grassmann algebra 13 - 1 : 



exp 



(Z\:Ait\t):\Z')=A{t,Z') 



x exp 



(64) 

(65) 
(66) 



vs. 



where £ are Grassmann vectors and |£) are Fermi coher- 
ent states. 

In order to use the above identities, it is necessary to 
transform e lc hc into a normal ordered form. Because 
h is Hermitian, it can be diagonalized by the unitary 
transformation such that h = UDU^ , where U is a uni- 
tary matrix and D is a diagonal one. With the canonical 
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{iDiY 



n=l 



= J] [l + (e iDi - 1)717.] 

= [] = e^'" 1 ^ :=: e *V"-J)e :, (67) 



Similarly, the projection onto the state which has N a 
electrons with spin a is defined as 

P"AX*) = ^ fj dcj){N„\f^{<P)\N^f^(<P) 

= ^ J J dfaxM, NJe"^*** (72) 
where g^f <T (4>,N a ) — e~ % ^ N " and 



we obtain 

gic^hc . Be . 

=: e &1 (e ih -i)c .. e c<Bc 



used. Then for any matrix £?, 










(73) 




I a = ^cr.T 


I 





o" 

/ 




(74) 



J ©tf*?2>y|€><€l : eW-W : |»|)<»/| : e* +B * : | 7 )< 7 | 

x exp[f(e ih - J)t/] exp 
x exp[r7B 7 ] exp 
V£VeVj\t) exp 



^ ~ 2^ ~ \ m 



1 



1 



»77 - TjW _ 2 77 



(7l 



7 1 TC l ~ 



x exp [e{e 4h (S + 7) - I} 7 ] 
1 1 



x exp 



£7 - 2 ee - 2 77 



(7l 



= . e ct[e ih (B+I)-I]c . 



(68) 



where e = e lh rj and 

= J] dI * de x = II drj x e- lh e lh d Vx = Vq. (69) 

Thus Eq. (61) is proven by taking B = —[21 + (n T — 

i)- 1 ]. 

^.4J Examples of quantum-number projector 
4-2.1 Particle-number projector 

Since the GBMC method is a grand canonical ap- 
proach, particle-number projection is needed when treat- 
ing the canonical ensemble. The projector onto a state 
with a given particle-number TV is defined as 

W = f dcj>{N\f M (<P)\N)^f u {<P) 
27T Jo 

= ^ J J dfaxtt, JV)e i£t ^W £ , (70) 

where gx{(j>, N) = e"^ and 

TV(</») = e i0E » e « t6i , e ih " w = e^I. (71) 



^.,2.,§ Total-spin projector 

The SU(2) symmetry of the total spin is recov- 
ered by summing up all the Euler angles in the spin 
space. 11 ' 16 ' 17 ) Thus, projection onto a given total-spin 
state is defined as (here we restrict ourselves to the case 
of S z = 0) 



Ps{S) = 



2S+1 
_ 25+1 

x P s (cosp)e ieths ^ a)e e ieths y w& e l& ' hs ^^ & , 

(75) 

where Ps(cos/3) is the S th Legendre polynomial and 



du(S,0\t s (u)\S,0)ifs(u) 

/ da d(3sm[3 / 
Jo Jo Jo 



t s (e,) = e iaS 'e i < lSv e* S ' 



hs z (a) 



hsM 



[3 





il 



il 




(76) 
• (77) 



Here, S a , (a = x,y,z) corresponds to the total a- 
component of spin : 



er a 



(78) 



Since the total-spin projection involves triple integrals, 
the computational cost is rather high. However, if one 
takes iVf projection and iVj, projection before the total- 
spin projection, the integrations about Euler angles a 
and 7 can be done analytically. Thus, the total-spin pro- 
jection is reduced to 

Ps(S) = £ df3g s ([3,S)e^ h ^ & , (79) 

where gs (P,S) = sin (5 P$ (cos (5) . 

4-2.3 Lattice- symmetry projector 

When the Hamiltonian is invariant under certain geo- 
metrical transformations, such geometrical symmetry is 
recovered by summing up all the transformations. For in- 
stance, when treating square lattice systems, they have 
Ci V lattice symmetry (see Fig. 10). By assuming that 
the square lattice lies in the x-y plane, the 7r/2 rotations 
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Fig. 10. Each element of C 4v group. 



For a simpler case, C4-symmetry projector is obtained 
by restricting a to the subgroup of C 4t) , i.e., 



9cM) T cA a ), 



{E, C 4 , C 2 , C 4 }, 



(90) 
(91) 



where the definition of gc 4 (a) is same as gc iv (a) within 
a = {E,C i ,C 2 ,Cl}. 



around the z-axis is achieved by the ^-component of the 
angular momentum L z . Let Tc 4 — e l % L * such that 



(80) 
(81) 



where C 4 (i) denotes a 7r/2 rotation around z-axis. From 
the above, e Jcthc 4 c is represented as 



3 ic t h. C4 c 



0'<t'),(»<t) 



aa'0jC 4 (i)- 



(82) 



Similarly, the x-y mirror transformation T„ d = e tcih "d c 
is defined as 



where <Td(i) denotes a x-y mirror transformation of i-ih 
site, i.e., (j d {x,y) — (y,x). Then, the representation of 
a d becomes 

~e i&]h °* & ] = <W<W (84) 

Although C iv group has other elements as in Fig. 10, 
all can be generated by C 4 and ad (see Table. II). For 
example, a x — C 4 ■ ad, i.e., 



C4 ■ ad 



= C 4 



-1 

1 









7. 


)=a x ( 









(85) 



Thus, the projection onto the C 4t) symmetry sector reads 



(86) 



where a = {E, C 4 , C 2 , Cf , a x ,a y , a d , a d >} and 
for s wave 

fcW = l(Va), 



for d x 2_ y 2 wave 



yG4 " v ; \ -1 (a = C 4 ,Ci,a d ,a d >) 
for d xy wave 

, x _ J 1 (a = E,C 2 ,a d ,a d <) 

j I -1 (a = C 4 ,Cf,<7 x ,a„) 



(87) 
(88) 

(89) 



Table II. The multiplication table of C 4v 
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c 4 


C 2 


cl 


&X 


Oy 


o d 


o d > 
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E 




c 2 


cl 


Ox 
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O d 
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C* 4 


C4 


c 2 


cl 
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°d 


Ox 


Oy 


C 2 


c 2 


cl 
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C 4 


Oy 
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o d > 


Od 


cl 


cl 
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C* 4 


C 2 


Od 


O d > 


Oy 


Ox 


Ox 


Ox 


°d 


(Jy 


°d> 
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c 2 


C* 4 


cl 


Oy 


o y 


o d > 


Ox 


<?d 


c 2 


E 


cl 


c 4 


°~d 


a& 


<Jy 


O d ' 


Ox 


cl 


c 4 


E 


c 2 


°d' 


o d > 




o-d 


(Jy 


c 4 


cl 


C 2 


E 



4-2.4 Total-momentum projector 

When the Hamiltonian has the translational symme- 
try, the total momentum must be conserved. A transla- 
tion by a lattice vector R is achieved by the operator 
T(R) = e lR ' ktot . Here, fe to t is defined by the Fourier 
transformation of the creation and the annihilation op- 
erators as 



k 



tot 



£*4a*> cL = ^E e ^4> ( 92 ) 



where i denotes a vector to the i-th site and N is 
the number of sites. The sum over k goes over all the 
points in the first Brillouin zone. The projection onto the 
Hilbert space with the total momentum K then reads: 

P K (K ) = ±Y,(K \f K (R)\K )lf K (R) 



R 



(93) 



R 



where gjc(R,Ko) = e lRK ° and the sum over R goes 
over all the lattice sites and hfc(R) is 



(h K (R)) ij = ±R.J2ke ik < i -»- 



(94) 



Then the matrix representation of e lcjfl 'c(R) c re ads 



(ia),(jo') 



= S a _ 

k 

= 5trcr' $i+R, j ■ 



(95) 
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4-3 Projected expectation values of observables 

From the previous subsection, one can define a general 
quantum-number projector as 



J dxg(x)f(x), 



(96) 



where T(x) is unitary and thus P* = P. For simplicity, 
we assume that the physical observable O commutes with 



P, i.e., 



P,0 



= . Then, the projected expectation 



value of the observable O becomes 



(O)p 



Tr 


PpPO 


Tr 


PpO 


Tr 


PpP 


Tr 


Pp 



(97) 



Here we use the projection property P 2 = P. Replac- 
ing the density-matrix operator by the sum over all the 
walkers p = J d\P(\)A(\) ~ A(AJ yields 



(0)p 



EjdxTr 


f(x)A(\JO 







_ J2iJ dxg(x)Q,i(x)Q(ni(x)) 
T,i I dxg(x)Qi(x) 

where f(x)A(\) = A(\(x)). 



(98) 



4-4 Post-projected sampling method 

In this section, we demonstrate the results of the 
GBMC-PS method. As we see in §3.4, the original GBMC 
method fails in reproducing the ground state of 4 x 1 lat- 
tice Hubbard model at U/t = 4 and n — 1. Here, we 
assess the accuracy obtained from the quantum-number 
projection to the density matrix obtained by the origi- 
nal GBMC method by following the idea of Assaad et 
a/. 16 - 17 ) Fi gure 11 shows the total energy compared with 
the exact diagonalization result. Here, we have projected 



the total z-component of spin S z — and the total spin 
S = 0. In this case, the exact ground state energy is 
-2.102748 whereas our data is -2.1037 ± 0.0010. The 
error is obtained by averaging the data over the imagi- 
nary time after the convergence. This result shows that 
the quantum-number projection method well reproduces 
the ground state energy, which is not obtained in the 
framework of the original GBMC method. This result is 
consistent with the result of Assaad et al. 16 ' 17 * 1 

However, the GBMC-PS method suffers from a slow 
convergence when the interaction strength U/t becomes 
larger. As is illustrated in Fig. 12, the energy in the case 
of U/t — 10 and n = 1 on 4 x 4 lattice under the full 
periodic boundary condition obtained by the GBMC-PS 
method is not yet converged with the ground state at 
t = 6. Since the strong on-site repulsion U/t prevents 
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Fig. 12. Energy as function of r at U/t = 10 and n = 1 on the 
4x4 lattice under the periodic boundary condition. Dashed line 
is the exact ground state energy obtained from the exact diago- 
nalization. Here, for a fast convergence, we use Green's function 
obtained by the Hartree-Fock calculation as an alternative start- 
ing point. This scheme for the fast convergence will be discussed 
in 84.8 in detail. 
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Fig. 11. Energy of the 4 X 1-site Hubbard model with U/t = 4 
and n = 1. Here, K = 0, S z = and 5 = projections are 
performed for r > 5. The integrals of S z and S are evaluated by 
the Gauss-Legendre quadrature 23 ) with the summation over 10 
points mesh each. 



the state from updating, the efficiency of the importance 
sampling of the original GBMC becomes worse with the 
increase of the interaction strength U/t. Thus, in the 
framework of the GBMC-PS method, it is difficult to 
store the samples which has a large overlap with the 
ground state. In the next subsection, we introduce the 
PR-GBMC method to overcome this slow convergence 
in the strong interaction regions. 

4-5 Pre-projected sampling method 

In this subsection, we introduce an alternative method 
for performing the quantum-number projection which 
is based on the importance sampling in combination 
with the quantum-number projection. This allows us to 
perform the sampling with the projected weight, which 
is more efficient than performing the sampling with 
the original weight. We call this pre-projection method 
the PR-GBMC method in contrast with the GBMC-PS 
method introduced by Assaad et al. 16 ' 17 ^ 



onto the state which has the total momentum K = 0, 
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For the pre-projected sampling, we rewrite Eq. (98) as 
(O)p 



where Q, = J dxg(x)Cli(x). Estimating (O) p is now 
reduced to the calculation of the weighted average of 
J dxg(x)Qi(x)0(ni(x))/tli with respect to the pro- 
jected weight Oj. 

The projected weight Sl^ stems from the projected 
density-matrix operator, i.e., 

Tr[P/3Pt] = j dxg{x)Tr [f (x)A(Sli, rn) 

= £ I dxgix^x) = ]T Q h (100) 

where f(x) e ic^h(x)c and 

fli det [(e lh ( x ' — I)nJ + /] . If the original samples 

have no overlap with the projected sector, 

Tr[PpPt] = & 



Oj / dxg(x) det (e 



ih(x) 



(101) 



becomes zero. However, from Eq. (47) the unprojected 
weight fli is always positive. It suggests that the fac- 
tor / dxg{x) det [{e 1 ^ - I)nf + l] in Eq. (101) causes 
the reduction of £li to tti when the original samples have 
small overlap with the projected sector. Empirically, we 
find that this reduction is realized by the cancellation of 
positive tti and negative flj, thus 12, is not positive def- 
inite. This is the source of the negative sign problem in 
the PR-GBMC method. In this case, we introduce the 
sign variable Si = ±1 by Cli = Si\Cli\, and the impor- 
tance sampling is performed with the absolute value of 
Oj, namely we calculate 

£ \q.\s^ dX9 ^ Qi ^ ^ ni ^ 



(102) 



Appearance of the negative sign is evaluated by calculat- 
ing the expectation value of the sign defined as 



(S) 



(103) 



Below we show results of the PR-GBMC method with 
the quantum-number projection K = 0, S z = and 
5 = 0. Figure 13 shows the energy of 4 x 1 lattice Hub- 
bard model with U/t = 4 and n = 1 calculated by the 
PR-GBMC method. Averaging the data over imaginary 
time gives the energy — 2.1026±0.0002 in agreement with 
the exact diagonalization result —2.102748. During the 
simulation, the average sign (S) is kept unity, i.e., there 
is no negative sign. We will discuss later the negative sign 
problem in more detail. 

One of the main advantages in the PR-GBMC method 
is that this method allows us to analyze directly the 




Fig. 13. Energy of the 4 X 1-site lattice Hubbard model with 
U/t = 4 and n = 1. Here, K = 0,S Z =0 and 5 = projec- 
tions are performed. 



change in the probability distributions caused by the 
quantum-number projection: 

dP{Q,n) 



dr 



jCP(Q,n) 

projection 



= £P{n,h), (104) 



where the projected variables are denoted by tilde. The 
reweighted distribution P cannot be calculated in the 
framework of the GBMC-PS method, because in the 
GBMC-PS method, the quantum-number projection is 
performed by reweighting the importance sampling of 
the original GBMC method, i.e., 

E* I dxg(x)Sli(x)0(ni(x)) 



(0) P = 



E 4 / dxg(x)Vti{x) 



Lift 



f dxg(x)Qi (x)O(n, (a:)) 



Lift 



f dxg(x)Q,i{x) 



sampling by Q,i N mc 



1 v^v f dxg{x)i1i (x)O(rii(x)) 



E, 



f dxg(x)Q,i(x) 



(105) 

where N mcs is the number of Monte Carlo sam- 
ples. As is seen from Eq. (105), what one can ob- 
tain by the GBMC-PS method is the distribution 
of J dxg(x)0.i{x)O{rii{x))/Q,i and J dxg(x)fli(x)/£li. 
Thus, the distribution of (O) p itself can not be calcu- 
lated by GBMC-PS method. 

We now calculate the projected probability distribu- 
tion P{N) to compare with the unprojected distribution 
P(N) obtained by the original GBMC method. Figure 14 
shows the distri bution of Green 's function. The abscissa 
TV represents ^Jj2ija n \ia) (ja) ^ or tnc GBMC method 

and 
where 



hf. w - / \ for the PR-GBMC method, 



n 



J dxg(x)Q( 



x)n(x)/\tt\. 



(106) 



As is illustrated in Fig. 14, the projected distribu- 
tion P{N) decays exponentially. This means that the 
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Fig. 14. Distribution of Green's function. The solid curve shows 
the distribution obtained by the PR-GBMC and the dashed 
curve shows the distribution obtained by the GBMC which is 
the same data as in Fig. 8(b). The data are obtained for the 
4 X 1-site Hubbard model with U/t = 4 and n = 1 at r = 20 
with 3.84 X 10 6 Monte Carlo steps. 



quantum-number projection actually reduces the phase- 
space. In the PR-GBMC method, the importance sam- 
pling is performed with respect to this reduced phase- 
space. Therefore, the convergence to the ground state is 
faster than the GBMC-PS method which is based on the 
importance sampling with the original weight. 



is more efficient. Actually the PR-GBMC method offers 
a better convergence at larger U/t. This is because the 
original GBMC sampling fails in making samples which 
have enough overlap with the ground state at relatively 
large U/t. This possibly causes a serious minus sign prob- 
lem. Thus when treating large U/t systems (typically 
U/t>4), the PR-GBMC method has to be employed. 
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Fig. 15. Energy of the 4x4 lattice with U/t = 4 and n = 1 as 
functions of r. The squares represent the result of the GBMC-PS 
method and the triangles represent the result of the PR-GBMC 
method. 



4.6 Comparison between GBMC-PS and PR-GBMC 

In the previous subsection, we have confirmed that 
the quantum-number projection changes the probability 
distribution and removes errors arising in the original 
GBMC procedure. Here, we make a comparison between 
the GBMC-PS and PR-GBMC methods to discuss their 
merits and demerits. 

Figure 15 shows the energy of the 4 x 4-site Hub- 
bard model with U/t — 4 and n — 1 under the pe- 
riodic boundary condition. As is seen from Fig. 15, 
the energy obtained by the PR-GBMC method con- 
verges with the ground state faster than that obtained 
by the GBMC-PS method. This comes from the differ- 
ence of the sampling procedure. In GBMC-PS method 
importance sampling is performed with respect to the 
unprojected weight fl = eT fa H(n)dr Thus the sam- 
pling depends only on the energy. On the other hand, 
PR-GBMC method makes use of the projected weight 
O = Q J dxg(x) det[(e ih - I)n T + J], which reflects 
not only the energy but also the overlap with the pro- 
jected sector. Therefore, the PR-GBMC method allows 
the convergence to the ground state at smaller r than the 
GBMC-PS method. However, in the PR-GBMC method 
we have to perform the projection for every sample, while 
the GBMC-PS method requires the projection only for 
accepted samples, which makes the computation time 
shorter. Empirically, we find that the energy obtained 
by both the GBMC-PS and the PR-GBMC methods con- 
verges with the ground state when the on-site interaction 
U/t is not too large, whereas the PR-GBMC method 



4-7 Negative sign problem 

Since the PR-GBMC method is based on samplings 
with respect to the projected weights fl, there appears 
minus sign if a sample has small overlap with the pro- 
jected sector. Figure 16 shows PR-GBMC results of 
the total energy and the expectation value of the sign 
(S) = |Oi|Si/Ei on the 4 x 4 lattice at U/t = 4 
and n = 1. As we see in Fig. 16, the average sign de- 
creases first, but it recovers along with the convergence 
of the energy. This means that if the overlap with the 
projected sector is small, the expectation value of the 
sign becomes small, while (S) recovers when the samples 
gain a large overlap with the quantum-number-projected 
state. We note that the r dependence of (S) is completely 
different from that in the conventional AFQMC method 
as well as from that in other methods, where (S) expo- 
nentially decreases to zero with increasing r. 

4-8 Acceleration of convergence 

In the previous subsections, we have proposed and 
studied the PR-GBMC method and examined the neg- 
ative sign problem. In the PR-GBMC method, the neg- 
ative sign appears when the overlap with the projected 
sector is small. To avoid this problem, instead of em- 
ploying no = \l as the initial condition, it is better to 
use Green's function obtained by the Hartrcc-Fock calcu- 
lation as an alternative starting point. Since the initial 
state is already that of the Hartree-Fock solution, the 
overlap with the ground state is expected to be relatively 
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of the projected distributions becomes slower with the in- 
crease of U/t. In this subsection, by comparing the U/t 
dependence of the energy convergence with that of the 
projected distribution, we discuss the applicable range of 
the PR-GBMC method. 

Figure 18 shows the energy of 4 x 4 lattice under the 
periodic boundary condition with n = 1 at U/t = 4, 10 
and 15. All the results are obtained by 1.6 x 10 6 Monte 
Carlo steps with the pre-projection at K = 0, S z = 
and S = 0. We employ the data at r = from the so- 
lutions of the Hartree-Fock calculation with U nt = U/2. 
Measurements are divided into 5 bins and the error bars 
are estimated by the variance among the 5 bins. As is 



Fig. 16. Energy and average sign of PR-GBMC results as func- 
tions of t. The squares represent the energy of 4x4 lattices under 
the periodic boundary condition with U/t = 4 and n = 1 and 
the triangles represent the average sign. The initial condition at 
t = is Green's functions for the infinite temperature. 



larger than n = \l- Therefore, the convergence to the 
ground state becomes faster and the expectation value 
of the average sign becomes stable. Figure 17 shows the 
total energy of 4 x 4 lattices under the periodic bound- 
ary condition with U/t = 4 and n — 1 together with 
its average sign. As is seen from Fig. 17, the total en- 
ergy converges already at r = 0.5 and the average sign 
is nearly unity in the whole range of r. 



-10 

-11 



-12 



-13 



-14 



it- 



<s> 

^ il4 ,' 1 '' 1 M 1-00 



Exact 

PR-GBMC(HF) ■ 

Sign A 

PR-GBMC • 



Hartree-Fock 
solution 



2 

x 



0.98 
0.96 
0.94 

0.92 
0.90 



Fig. 17. Energy and average sign of PR-GBMC results as func- 
tions of t. The squares represent the energy of 4 X 4 lattices 
under the periodic boundary condition with U/t = 4 and n = 1 
started from the Hartree-Fock solution (open star) and the trian- 
gles represent the average sign. For comparison, the PR-GBMC 
result obtained from the infinite temperature at r = is plotted 
(circles). 



4.9 Applicability of PR-GBMC method 
4-9.1 U/t dependence of convergence 

In the PR-GBMC method, the distributions of the 
phase-space variables P(£l) and P(N) are transformed to 
the projected distributions P(Cl) and P(N) which decay 
faster than the original distributions. However, the decay 
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Fig. 18. Energy as functions of t at (a) U/t = 4, (b) U/t = 10 
and (c) U/t = 15 for the 4 X 4-sitc Hubbard model under the 
periodic boundary condition . Dashed line is the exact ground 
state energy obtained from the exact diagonalization. 
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seen from Fig. 18, the convergence with the ground state 
becomes slower with the increase of U/t and at U/t = 15, 
the energy obtained by the PR-GBMC method does not 
yet converge with the ground state at r = 2 because 
of the slow convergence. At U/t = 15, since the energy 
seems to show further decrease beyond r = 2, further 
time evolution is required whereas the statistical error 
becomes large. This large statistical error comes from 
the fast diffusion of the probability distribution of the 
sampling weight f2. Since the definition of the weight is 
given by fi = e~ & H ( n ') dT ( the weight O grows expo- 
nentially with r. Thus, the phase-space to be sampled 
becomes larger if the convergence with the ground state 
becomes slower because of the strong interaction. This 
requires much more computation time, which determines 
the practical limitation of the PR-GBMC method. 

To confirm this, we calculate the distributions of 
the projected phase-space variables f2 and N — 

Figures 19 and 20 show the inte- 



Q(Q) 



grated distributions Q(f2) and Q(N) defined by 

Q(fi) = 1-1 P(Q')dfl', 
Jo 

N 



Q{N) = 1 - / P(N')dN' 



(107) 



(108) 



As we see in Fig. 19, the Monte Carlo step dependence 
of Q(Cl) implies the existence of the cutoffs in the distri- 
butions of fl for all U/t. Thus, the boundary terms with 
respect to the projected weight 57 do not appear to ex- 
ist for all U/t. However, the phase-space of the projected 
weight Cl becomes larger with U/t and at U/t = 15, and a 
distinct plateau structure with steps caused by the lack 
of the large Ct samples is seen in the tail of the distri- 
bution (see Fig. 19. c). This means that in the case of 
U/t = 15, 1.6 x 10 6 Monte Carlo steps are not enough 
for the accurate sampling of events at large weight CI. 
Since the samples with large weight seldom appear, but 
contribute to lowering the energy, these samples cause 
the spike structure in the distribution of the energy (see 
Fig. 21). Thus, the statistical error of the energy becomes 
larger. 

The slow convergence of the distribution due to the 
strong interaction is also observed in the distribution of 
the projected Green's function Q(N). As is seen from 
Fig. 20, the convergence of the distribution tail becomes 
slower as U/t increases. Although the energy converges 
with the ground state at U/t = 10 as in Fig. 18. b, the 
plateau structure in Q(N) visible at U/t = 10 (Fig. 20. b) 
signals the slow convergence in the PR-GBMC method. 
However, the convergence of the energy at U/t = 10 
shows that a small plateau structure arising in the tail 
part of Q(N) does not yet cause a bad effect on the 
convergence of the energy. At U/t = 15, the distribution 
tail of Q(N) shows no indication of the existence of the 
cutoff at least up to 1.6 x 10 6 Monte Carlo steps, which 
indicates that the number of Monte Carlo steps 1.6 x 10 6 
is not enough at U/t = 15. 

These slow convergences of the distributions at large 
U /t come from the expansion of the phase-space to be 
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Fig. 19. (color online) Integrated distribution of projected weight 
Q at r = 2 for the 4 X 4-site square lattice under the periodic 
boundary condition with (a) U/t = 4, (b) U/t = 10 and (c) 
U/t = 15. In all the panels, red (light), green (long dashed), 
blue (dashed), pink (dotted) and black (dark) curves represent 
the distribution obtained by 3.20 X 10 5 , 6.40 X 10 5 , 9.60 X 10 5 , 
12.8 X 10 5 and 16.0 X 10 5 Monte Carlo steps, respectively. 



sampled. From Eqs. (38-40), the drift term and the diffu- 
sion terms in Langevin equation (44) are proportional to 
U and y/U, respectively. Thus the diffusion speed of the 
distributions becomes faster with the increase of U/t, 
which means that the phase-space to be sampled be- 
comes larger. This larger sampling space requires larger 
Monte Carlo steps and causes an insufficient sampling 
at large U/t. Therefore, it is advisable to monitor the 
convergence of the distributions Q(Cl) and Q(N) to en- 
sure that the number of Monte Carlo steps is sufficiently 
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Fig. 20. (color online): Integrated distribution of projected 
Green's function N at r = 2 for the 4 X 4-site square lattice under 
the periodic boundary condition with (a) U/t = 4, (b) U/t = 10 
and (c) U/t = 15. In all the panels, red (light), green (long 
dashed), blue (dashed), pink (dotted) and black (dark) curves 
represent the distributions obtained from 3.20 X 10 5 , 6.40 X 10 s , 
9.60 X 10 5 , 12.8 X 10 5 and 16.0 X 10 5 Monte Carlo steps, respec- 
tively. 



large. 

4-9.2 Size dependence of convergence 

In this subsection, we analyze the size dependence of 
the convergence by employing the PR-GBMC method 
with K = 0, S z = and 5 = projections. Figure 
22 shows the energy of 6 x 6, 8 x 8 and 10 x 10 lat- 
tices at U/t = 4 and n — 1. As is seen from Fig. 22, 




Fig. 21. Distribution of energy at U/t = 15 for 4 X 4-site Hubbard 
model under the periodic boundary condition at n = 1. The data 
are obtained from 1.6 X 10 6 Monte Carlo steps at r = 2 



all the simulation results converge with the ground state 
energy obtained by the AFQMC method with the Trot- 
ter discretization of At = 0.025. As we see in Fig. 23, in 
accordance with the convergence of the energy, the inte- 
grated distributions of the projected weight Q((l) at any 
size show no distinct signal of insufficiency in sampling 
which is observed in the tail part of Q(Cl) at U/t = 15 
on 4 x 4 lattice as illustrated in Fig. 19. c. The existence 
of the cutoffs in the integrated distributions of Green's 
function Q(N) also supports the fact that the number of 
Monte Carlo steps 2.56 x 10 4 is sufficient in these systems 
(see Fig. 24). 

Here, to confirm the convergence, we have calculated 
not only the ground state energy but also several physical 
quantities. First, we evaluate the equal-time spin and 
charge correlations defined by 



S(k) = ^J2 elk ' {r '~ r3) &-S- 



31 ' 



(109) 



N(k) 



£ J2 e*^-^((n 4T + ftii)(n iT + h n )). 



(110) 



As we see in Fig. 25, for all the system sizes, the peak 
values of the spin correlation S(tt, it) obtained by the PR- 
GBMC method are consistent with the values obtained 
by the AFQMC method. The peak values of S(n, ir) as 
well as the charge correlation N(ir, ir) are listed in Table 
III. A small discrepancies observed in N(tt, n) may be 
attributed either to the effect of finite intervals in the 
imaginary time step of AFQMC or to the slight insuf- 
ficiency of r taken for the ground-state average in the 
PR-GBMC calculation. 

Next, we show the superconducting correlation defined 

by 



(111) 



where P a (r) denotes the equal-time pairing correlation 
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Fig. 22. Energy as functions of r at U/t = 4 and n = 1 for (a) 
6x6, (b) 8x8 and (c) 10 X 10 lattices under the periodic boundary 
condition. Dashed lines are the ground state energy obtained by 
the AFQMC method. Error bars are as large as the symbol size. 
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defined as 

1 N 

P ° W = 2N E< A ^) A «( l + r ) + (i + 0), 

i=l 

(112) 



where r is the distance from the i-th site and A a is the 
superconducting order parameter. The latter is defined 
as 

A a («) = -j=^fa{r){Ci-\Ci+rl - c a c 4+rT ), (113) 

where f a (r) is the form factor of the pairing correlation 
defined as 



fi s (r) = A5 rxfi 8 ryfl , 

hs{r) = 5 r o{S Tmj \ + 5r a ,-l) 



(114) 



Fig. 23. (color online): Integrated distribution of projected 
weight Q obtained att//t = 4n = l and r = 1 on (a) 6x6, (b) 
8x8 and (c) 10 X 10 lattices under the periodic boundary con- 
dition. In all the panels, red (light), green (long dashed), blue 
(dashed), pink (dotted) and black (dark) curves represent the 
distribution obtained by 0.512 X 10 4 , 1.024 X 10 4 , 1.536 X 10 4 , 
2.048 X 10 4 and 2.560 X 10 4 Monte Carlo steps, respectively. 



+ *r.,o(*r 1) ,l + *r s ,-l). (H5) 
fsd(r) = £r B ,u(<W + S Tw ,-l) 

- S rx ,o(Sr vl l + S ry ,-i), (116) 

/3s W = S ri ,l(Sr y S + $r y ,-l) 

+ 5r.,-l(*r„-l + 5p v ,l). ( 117 ) 
/3d M = S riil (S ry:1 - S ry -x) 
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Fig. 24. (color online): Integrated distribution of Green's func- 
tion N obtained at U/t = 4 and r = 1 on (a) 6 x 6, (b) 8 x 8 
and (c) 10 X 10 lattices under the periodic boundary condition. 
In all the panels, red (light), green (long dashed), blue (dashed), 
pink (dotted) and black (dark) curves represent the distribution 
obtained by 0.512 X 10 4 , 1.024 X 10 4 , 1.536 X 10 4 , 2.048 X 10 4 
and 2.560 X 10 4 Monte Carlo steps, respectively. 



+ S rie ,-i(5 r -x - 5 r i), 



(118) 



where Sij is Cronecker's delta. The suffices a — 
Is, 2s, 2d, 3s and 3d represent the on-site s-wave, the 
extended s-wave directing along the x and y axes, the 
d x 2_ y 2-wave, the extended s-wave along the diagonals 
and the d xy -wave, respectively. Numerical results for all 
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Fig. 25. Peak values of spin correlation S(n, tt) on (a) 6x6, (b) 8 X 
8 and (c) 10 X 10 lattices under the periodic boundary condition 
at U/t = 4 and n = 1. In all the panels, solid and dashed lines 
represent the results obtained by the AFQMC method and their 
error ranges, respectively. 



the quantities defined above are listed in Table III. Here, 
for comparison with the other method, we demonstrate 
the numerical results obtained by the AFQMC method 
with the Trotter discretization of At = 0.025 and the 
number of Monte Carlo steps N mcs = 10 4 . 

We have also calculated the momentum distribution 
defined by 

In Table IV, we show the numerical results for the mo- 



20 J. Phys. Soc. Jpn. 



Full Paper 



Online-Journal Subcommittee 



Table III. Comparison between PR-GBMC and AFQMC meth- 
ods. For all the parameter sets, we employ the quantum-number 
projections onto the total momentum K = 0, total z-componcnt 
of the spin S z = and the total spin S = 0. All the simulation 
results are obtained by 2.56 X 10 4 Monte Carlo steps under the 
periodic boundary condition. 



6x6,!//i = 4,n = l 


PR-GBMC 


AFQMC 


Energy 


-30.87 ± 0.02 


-30.87 ±0.03 


S(n,n) 


1.43 ± 0.01 


1.43 ±0.03 


N(tt,it) 


0.409 ± 0.001 


0.403 ± 0.001 


Sis 


1.636 ± 0.005 


1.610 ±0.004 


S2s 


1.188 ± 0.004 


1.184 ±0.001 


S2d 


1.101 ± 0.007 


1.116 ±0.018 


S 3 s 


0.539 ± 0.004 


0.532 ±0.003 


S3s 


0.397 ± 0.003 


0.399 ± 0.012 


8 x 8,U/t = 4,n = 1 


PR-GBMC 


AFQMC 


Energy 


-55.01 ±0.03 


-55.09 ±0.06 


S(tT, 7r) 


2.08 ± 0.04 


2.05 ±0.04 


jV(7r,7r) 


0.420 ± 0.003 


0.412 ± 0.002 


Sis 


1.681 ± 0.011 


1.650 ±0.001 


S2s 


1.191 ± 0.001 


1.190 ±0.001 


S2d 


1.077 ±0.004 


1.110 ±0.018 


Szs 


0.515 ± 0.004 


0.505 ±0.002 


Sis 


0.457 ± 0.005 


0.457 ±0.007 


10 x W,U/t = 4,n = 1 


PR-GBMC 


AFQMC 


Energy 


-86.25 ±0.10 


-86.11 ±0.03 


S(tT, 7r) 


2.88 ± 0.03 


2.88 ±0.06 


jV(7T,7r) 


0.424 ± 0.002 


0.417 ±0.001 


Sis 


1.694 ±0.006 


1.668 ±0.004 


S2s 


1.195 ± 0.005 


1.192 ±0.001 


S2d 


1.076 ± 0.005 


1.111 ±0.018 


Szs 


0.510 ± 0.003 


0.507 ±0.002 


S\is 


0.471 ± 0.005 


0.459 ± 0.027 



Table IV. Momentum distribution n(k x ,k y ) obtained by PR- 
GBMC and AFQMC method. The results have been obtained 
from the same simulation as those of Table HI 



6 x 6, U/t = 4, n = 1 


(k x , k y ) 


PR-GBMC 


AFQMC 


(0,0) 
(w/3,w/3) 
(2tt/3, 2tt/3) 

(-7T, 71") 


0.967 ± 0.001 
0.928 ± 0.001 
0.072 ± 0.001 
0.033 ± 0.001 


0.966 ±0.001 
0.926 ± 0.001 
0.074 ± 0.001 
0.034 ± 0.001 


8 x 8, U/t = 4, n = 1 




PR-GBMC 


AFQMC 


(0,0) 
(w/4, W /4) 
(w/2,w/2) 
(3tt/4, 3tt/4) 

(-7T, 71") 


0.967 ±0.001 
0.951 ± 0.001 
0.500 ± 0.000 
0.049 ± 0.001 
0.033 ± 0.001 


0.966 ±0.001 
0.950 ± 0.001 
0.500 ± 0.000 
0.050 ±0.001 
0.034 ± 0.001 


10 x 10, U/t = 4, n = 1 


(h X , ky ) 


PR-GBMC 


AFQMC 


(0,0) 

(7r/5,7r/5) 

(2tt/5, 2tt/5) 
(3tt/5, 37I-/5) 
(4tt/5, 4tt/5) 
(n,n) 


0.968 ± 0.001 
0.959 ± 0.001 
0.883 ± 0.003 
0.117 ±0.003 
0.041 ± 0.001 
0.032 ± 0.001 


0.966 ±0.001 
0.957 ±0.001 
0.878 ± 0.001 
0.121 ±0.001 
0.043 ± 0.001 
0.033 ± 0.001 



mentum distribution along the line from (0, 0) to (tt } it) 
in the momentum space. 

In the numerical results shown in Tables III and IV, 
all the numerical data are essentially consistent each 
other. We also show a hard test of the convergence for 
the square-lattice Hubbard model with the next near- 



est neighbor transfer t' — 0.5i in one of the diagonal 
direction. This is nothing but the anisotropic trianglar 
lattice. This geometrically frustrated lattice structure 
generates a serious difficulty in various simulations. The 
PIRG method offers the only available results. 8 ' 9 ^ Here 
PR-GBMC results for the 6 x 6 lattice with the peri- 
odic boundary condition at U/t = 4 shows the ground 
state enrgy E — —32.54 ± 0.05 and the double occu- 
pancy (D) = 0.164 ± 0.001, which is favorably com- 
pared with the PIRG result of E = -32.6 ± 0.1 and 
(D) = 0.168 ± 0.006. This parameter value is, as is es- 
timated from the systematic studies of the PIRG result, 
just near the Mott transition point and provides us with 
a severe numerical benchmark. From the analysis of the 
size dependence, we conclude that the convergence and 
the applicability of the PR-GBMC method is not re- 
stricted by the system size. 

5. Summary and Discussion 

In this paper, we have reexamined the Gaussian-basis 
Monte Carlo method (GBMC) proposed by Corney and 
Drummond. 12, 13 - ) This method does not suffer from the 
minus sign problem for any Hamiltonian composed of up 
to two-body interactions (see Appendix A and B). How- 
ever, the original method often shows systematic errors 
especially in the low-temperature region. We have elu- 
cidated how the systematic errors come from the slow 
relaxation caused by the trap in the excited states. 

To overcome the systematic error, we have improved 
the quantum-number projection scheme proposed by As- 
saad et al. 16 ' 17 ) to make it possible to combine the projec- 
tion procedure in conjunction with the importance sam- 
pling of the original GBMC method. This method allows 
us to project out the excited states and improve the be- 
havior of the probability distributions, which makes it 
possible to widen the region of tractable parameters in 
the GBMC method. 

We have also discussed the applicability of our algo- 
rithm. In the PR-GBMC method, the convergence with 
the ground state becomes slower with the increase of the 
interaction strength U/t. This slow convergence is caused 
by the increasing inefficiency in the importance sampling 
procedure owing to the barrier in the phase-space com- 
ing from the strong on-site repulsion. Nevertheless, good 
convergence at U/t = 10 indicates a better efficiency 
of the PR-GBMC method compared to other numerical 
methods such as the AFQMC method. The system size 
dependence up to 10 x 10 lattice shows no distinct symp- 
tom of the slow convergence in physical quantities with 
the increase of the system size. Thus the slow conver- 
gence does not restrict the applicability of this method 
for large system sizes at least up to 10 x 10 lattices. 

In the large U/t region, the slow convergence is tightly 
associated with the exponential broadening of the dis- 
tribution of the weight Q, which causes an expansion of 
the phase-space to be sampled and results in inefficient 
importance sampling. Especially, the exponential growth 
of f2 and corresponding expansion of the phase-space of- 
ten results in the lack of rare event with large fi. The 
lack of Monte Carlo samples with large causes a spike 
structure in the tail part of the distribution P(f2) and 
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causes a plateau structure in the tail part of the inte- 
grated distribution Q(Ct). Despite samples with large fl 
contribute to physical quantities, such samples seldom 
appear at large U/t, which causes large statistical er- 
rors. Thus, with the increase of U/t, the convergence of 
the energy becomes slower accompanied by the increase 
of the statistical errors. If that is the case, the compu- 
tation time required for the convergence of the energy 
with required statistical errors goes beyond allowed com- 
putation time, which determines the practical limitation 
of the PR-GBMC method. Therefore, it is advisable to 
monitor the convergence of the distribution to evaluate 
the number of Monte Carlo steps needed for the conver- 
gence, especially when one calculates large U/t regions. 

With the inspection of the large U/t systems and the 
system size dependence, the PR-GBMC method offers 
a powerful tool which can be applied to the systems in 
several cases beyond tractable parameters of the conven- 
tional numerical methods such as the AFQMC and PIRG 
methods and is at least complementary to the existing 
methods. 
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Appendix A: Langevin Equations for a General 
Hamiltonian 

In this Appendix, we show how to construct a Gaus- 
sian representation for a general Hamiltonian. Here for 
simplicity, we treat a general number-conserving Hamil- 
tonian given by 

ija ijklaa' 



ijklaa' 



(AT) 



Using the operator identities in Eqs.(29)-(31), one-body 
term Hij„ gives a contribution Af- a to the equation of £1 
and a contribution A™^ to the drift term of the Green's 
function: 



ija 



ija 



where x (and y) denotes the site and the spin, i.e, x 
(i,a). The concrete expressions of Af- a and A™^" are 



Ajir _ tij<J n (icr),{jo)-, 

A ^xy tlj(j r , , 

A iju = — 7T [nx,(j<r){(>(i a ),y - n{i*),yi 



(A-3) 



}n ] . (A-4) 



Next, to guarantee the positive diffusion, we change 
the two-body term Hijkiaa' by adding the identity 

ft (»a),OV)~^'"(^).0>) =0 as 

Hijkiaa' = V r feYa'"(»(7),0-<T)«(fc<7'),(i<7') 



2 



+ 



TO I 



(A-5) 

The two-body term then gives a contribution A^- kltJtT , to 
the equation of fi, a contribution A™/,", , to the drift 

1 y i J Kiaa 

term of n xy and the contributions B^iaa 1 , C\jkiaa' to 
the diffusion term of n T 



L xy • 



(A-6) 



ijklaa' 



dn xy = Yl [ A W*„idT + B^,dW%l laa , 



ijklaa' 



+C?AH,dW, 



(2) 



. (A-7) 



y ijklaa' ijklaa' 

where the Wiener increment dW satisfies 

( dW ijllaa' dW i'j'k'l' m ') = dTdrr'Su'SjfSkk'Sll'Sa^a'v' 

(A- 

The concrete expressions of each term become 



4 s 



n 

ijklaa' 



V kla' 



in 



(ia),(j<r) n (ka'),(la>) 



-"(*a),(J<7') rl (fea').C7>)} 
+n(ia),(la')5kjSaa' + n(ka'),(ja)&il&aa'~\ , (A-9) 

v lja , 

A 7jklaa> = Y~ { n x,(la')(S(ka'),y-n { ka'),v) 

+ {&x,(la') - n x .(la'))n(ka'),y} 
+ n (ka>),(la') { n x,(ja)(S(ia),y ~ n (ia),y) 
+ ( S x,(ja) - n x ,(ja))n(ia),y} 

- { n (ia),(la>) - ^{ia),{la')} 

x { n x,{ja){5(ka'),y - n(ka'),y) 

+ ($x,(ja) - n x,{ja)) n {ka'),y) 

- {n(ka'),(ja) - 2<W),0>)} 
X { n x,(la>){5(ia),y - n {ia).,y) 

+ {Sx,{la>) - "-a;.(/<j')) n (^)-y}] ' 

(A-10) 



and 



ijklaa' 



|U yCT ,| 

k o' [ n x,(ja){5(ia),y ~ n (ta),y} 



-sign(VX,)n x ,(la>){$(ka>),y ~ n(ka'),y} 



(All) 
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-sign(y fe Y^){5 Xi(;<T / ) - fij,^)}^^),,, 

(A-12) 

In all, the general Hamiltonian (A-l) gives the Langevin 
equations of Ito-type 



dfl = 



xy 



ijff "T" ^ijklaa' 



ijklaa' 



Qdr, 



dr 



(A-13) 



+ \^ ijklaa' ^ijklaa' + ^ ijklaa' ^ijklaa' 

(A-14) 



ijklaa' 



Appendix B: Completeness and Positivity of 
the Gaussian Basis 

In order to establish a phase-space representation 
based on the Gaussian operators, one must show that 
they form a complete basis for the class of density-matrix 
operators which we wish to represent. Furthermore, the 
expansion of any density-matrix operators in terms of the 
Gaussian operators must involve the positive coefficients. 
For the self-contained description of this paper, we prove 
this 'positive completeness' by following the idea of J. F. 
Corney. 24 ) 

To prove this 'positive completeness', we will relate 
the Gaussian operators to the number-state projection 
operators which form a complete basis set and we will 
show that any term which appears in a number-state 
expansion of a density-matrix operator can be written as 
a sum of Gaussian operators with positive coefficients. 
Here, we prove only the number-conserving case in this 
Appendix, a proof for the most general Gaussian can 
be similarly constructed. 24 ) To this end, we rewrite the 
Gaussian operator in terms of \jl which is defined as 

/x = 21 + (n T - I)- 1 . 



(B-l) 



A number-conserving Gaussian operator A(n) is then 
written as 

A = dct(/-n) :e - &t [ 2 ' r +(™ T -' r r 1 ] & : 



1 



dct(2J - fi) 



M 



=■■ n( i -^W^):/det(2J-/x). (B-2) 

ij = l 

To avoid singular behavior, the limit of any rijj — > 1 is 
taken only in the normalized form of the Gaussian. 

B.l Number-state expansion 

Let n be a Fermionic occupation number vector n = 
(ni,ri2, • • • , nu), where = 0, 1, then a complete set of 
Fermionic number-state is represented by {|n)}, where 



{\ri}} runs over all the 2 M permutations. This set de- 
fines a complete operator basis of dimension 2 2M , and it 
enables us to expand the density-matrix operator as 

p = ^^\n)(n\p\m) (m| 

n m 

= Pnrn\n) (m\ = ^Y Pnrn, (B-3) 

n rh h rh 

where we impose a number-conserving condition J^i n i — 
J2i m i- From the positive definiteness of the density- 
matrix operator, all the diagonal density-matrix elements 
are real and positive: pnn > 0. Here we require addition- 
ally that pfin < 1 for the normalization. 

Since a density-matrix operator is Hcrmitian, it can 
be always diagonalized. Let \*S? k ) — J^CWIn) are * ne 
eigenvectors and P k the corresponding positive eigenval- 
ues of the density matrix, we can write 

p = £i\|* fc ><* fc | =J2^2P k C kH CUn)(m\. (B-4) 

k k h,rh 

Thus the coefficients of the number-state expansion can 
be represented as 



P. 



hrn ^ ' Pk^kn^kfh ^ ' ^kn^krhi 



(B-5) 



where D k n = \fP~k~Ckh- By using a Cauchy-Schwartz in- 
equality, the magnitude of these coefficients is given by 

\Pnm\ 2 = Y, ffe.n-DfeA-Ofe'n-^fc'm 
k,k' 

< ^ D k,nD*kn D k'ihDk'm 
k,k' 

= PnnPmrn < 1- (B-6) 

Thus the magnitude of any off-diagonal element is 
bounded, Conversely, any diagonal element is at least 
as large as the squared magnitude of anything else on 
the same row or column: 



\Pnfh\ — PnhPrhrh _ Pnn 
|2 



(B-7) 

|Pmn| — PrhfhPhn — Pnn- (B'8) 

Thus we obtain the lower limit of the diagonal elements: 

Pnn ^ max (|pnm| , |Pmn| ) 



> 



2(M 



— Y 



\Pnrh | 



> 



2(2 M - 1) ^ 



( | Pnrh j H~ |/?mn| 



(B-9) 



The number-state expansion of the density-matrix oper- 
ator can then be written as 

1 



P = Y A Pnn\n)(n\ + YY 



f-', 2(2 M - 1) 

n rh^n 



(B-10) 



J. Phys. Soc. Jpn. 



Full Paper 



Online-Journal Subcommittee 23 



where 



Pnrn{p) = |n)(n| + \fh) (m\ + 



2(2 M - 1) 

p^~ 

' urn 



n){m\ (B-ll) 



Apnn — Pnn 2{2 M — 1) ^ ' (|PrJm| H~ |Pmn| ) ^ 0. 

(B-12) 



Since all the coefficients of the new expansion (B-10) are 
positive, it is sufficient to prove that each operator in the 
expansion (B-10) can be written as a Gaussian or as a 
positive sum over Gaussians. 

B.2 Diagonal number-state projector 

First we show that the diagonal number-state projec- 
tor \n)(n\ in the new expansion (B-10) with positive co- 
efficients Apfin corresponds to a Gaussian operator. For 
individual ladder operators, one has the well-known iden- 
tities: 



|0>i<l|i, 



|l>j<0|i- 



(B-13) 



If we set riij = rijSij in Eq. (B-l), then /i^ = (5^(1 — 
2rij)/(l — rij) and the Gaussian operator reduces to 



M 



A=: II (l-^A) :/det(2J- M ) 



M 



II [(l-ni)MJ+»j^ 



M 



= Hl(l-n J )\0) J (0\ j +n j \l) j (l\ j }. (B-14) 

i=i 

Thus, if nj is chosen as or as 1, the Gaussian oper- 
ator itself can be regarded as a diagonal number-state 
projector: 



\n)(n\ = A(nj5ij). 



(B-15) 



B.3 Off- diagonal number-state projector 

Second we show that the mixed projector Pf L m{p) m 
the expansion (B-10) corresponds to a positive sum over 
Gaussians. Consider a Gaussian operator with 



njSij 



E 

k=i 



v k 5 



i.r(k) Oj.s(k) i 



(B-16) 



where the rij for each j is each either or 1 and the loca- 
tions (r(k),s(k)) of the Nk nonzero off-diagonal elements 
Vk satisfy r(k) ^ r(k'), s(k) ^ s(fc') and r(k) ^ s(k') for 
any k,k'(k ^ k'). In other words, if there is a nonzero 
element in the off-diagonal location (i,j), then there will 
be no other element in the i-th. row and j-th column and 



none in the j-th row and i-th column, i.e., 



n = 





o 






o 









rh 





riy 


0| 

























Tlj 


0| 






Si 






_p_ 


n M 



(B-17) 



This structure means that det(J — n) = YljLii^- 



and 



1 - 2rij 

N k 



fe=i 



r(k)Oi yS (k) 



Vk 



(1 -n r(k ))(l -n s(k ))' 



(B-18) 



Here again, we take the limit nj — > 1 only in the normal- 
ized form of the Gaussian to avoid the singularity. With 
these conditions, the Gaussian operator reduces to 



A=: 1] !',>'>.) :/det(2J-/x) 

l 

[< 



i,J = l 
M 



n k 1 -*] 



rijUjbj 



1 + 



Vkbl {k) br{k) 



(1 



V(fc) 



)(l-»«(fc)) 



><n 

fe=i 

M 

=: J| [(1 -^OIO^^O^- +« J -|1) J -<1| 
-r-r v kbl (k) b r(k ) 



(B-19) 



where the sum over if is the sum over all the possible 
subsets of {1, 2, • • • , Nk} and has 2 Nk terms. For the j-th 
mode, if there is a A: such that s(k) = j, then the diagonal 
number projector for the j-th mode is replaced 

by the off-diagonal |l)j(0|j, or if r(k) = j, the conjugate 
projector |0)j(l|j is created. Thus the Gaussian operator 
can be represented by a sum over number-state projec- 
tors: 



A = £(±) n «* |<s 

K kGK 



MK) 
l s(K) 



= K(n), (B-20) 
where the j-th element of the vector ffS^ is defined as 

(B-21) 



1 j = r(k),3k e K 
j = s(fc),3fc e if 
Hi otherwise 



A minus sign appears if an odd number of transpositions 
are required to put all the annihilation and the creation 
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operators in a canonical order. In this sum over projec- 
tors, the diagonal projector \n)(n\ is contained with coef- 
ficient 1. In the sum, there also exists the projector that 
transposes all the 2Nk specified modes with coefficient 
Y\k=i v k- The sum also contains projectors that trans- 
pose only subsets of these modes. In total, there are 
terms. 

Next by adding other Gaussian operators, we eliminate 
all the intermediate terms from the sum in Eq. (B-20) 
and leave only K — N k terms. First we add the N k 
Gaussians with one fewer off-diagonal element in the n 
matrix, to cancel the projectors that transpose N k — 1 
modes. Second we add jv fe C2 Gaussians with two fewer 
off-diagonal elements, to cancel the projectors that trans- 
pose Nk — 2 modes. This process is repeated until all the 
Sfc^i" 1 N u Ck = 2 Nk — 2 intermediate terms are removed. 
Finally, we obtain 



^ K,(n) = (2 N «-l)\n)(n\ 



v'Cv 



±n^ 



-.s(K)\ /-.r(K) 



k=l 



L r(K)/ \ n s(K) 



(B-22) 



m = rn^^y respectively, we 



where the sum indexed by subsets of u refers to the sum 
described above. 

By adding K v i (n) and K v i (m) , with different diagonal 
components n — and 
obtain 

Pfim(p) 
where 



2(2 M - 1) 
+ |m)(m| H \n)(m\ 



P- - 
1 inn 



2 N » 



(B-23) 



= ± 



Pnrn (2 M -1)(2^-1) 



(B-24) 



Thus it is proven that Pn m {p) can be represented by a 
positive sum over Gaussians and hence it is shown that 
any number-conserving density-matrix operator can be 



expanded by the Gaussian operators with positive coef- 
ficients. 
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